Learn how to use medfate for forest water balance simulations.
We are given all the necessary data, bundled in a single list:
fb <- readRDS("../StudentRdata/fb_data.rds")
names(fb)
## [1] "siteData" "treeData" "shrubData" "customParams"
## [5] "measuredData" "meteoData" "miscData" "soilData"
## [9] "terrainData" "remarks" "sp_params" "forest_object1"
Element fb$treeData is already in the appropriate format
for forest objects. Therefore, we can easily plug the tree
data into a forest object:
fb_forest <- emptyforest()
fb_forest$treeData <- fb$treeData
and examine its characteristics:
summary(fb_forest, SpParamsMED)
## Tree BA (m2/ha): 24.4861797 adult trees: 23.8297536 saplings: 0.6564261
## Density (ind/ha) adult trees: 3360 saplings: 1248 shrubs (estimated): 0
## Cover (%) adult trees: 100 saplings: 27.9158707 shrubs: 0 herbs: 0
## LAI (m2/m2) total: 2.7 adult trees: 2.4418971 saplings: 0.2581029 shrubs: 0 herbs: 0
## Fuel loading (kg/m2) total: 0.8739882 adult trees: 0.7934566 saplings: 0.0805316 shrubs: 0 herbs: 0
## PAR ground (%): NA SWR ground (%): NA
Note, incidentally, that LAI estimates come from the sum of column
LAI, without any use of allometric relationships.
A data frame with soil physical attributes are defined in:
fb$soilData
## widths clay sand om bd rfc
## 1 300 39 26 6 1.45 50
## 2 700 39 26 3 1.45 65
## 3 1000 39 26 1 1.45 90
## 4 2500 39 26 1 1.45 95
We need, however, to initialize a soil object to
complete the hydraulic parameters: .code80[
fb_soil <- soil(fb$soilData)
From which we can estimate the water holding capacity and extractable
water capacity for each layer (in mm) using functions
soil_waterFC() and
soil_waterExtractable():
soil_waterFC(fb_soil)
## [1] 58.32111 92.92727 37.30063 46.62578
soil_waterExtractable(fb_soil)
## [1] 26.06443 41.96683 16.97066 21.21332
So the sums would be:
sum(soil_waterFC(fb_soil))
## [1] 235.1748
sum(soil_waterExtractable(fb_soil))
## [1] 106.2152
The same information can also be found in the output of
summary().
summary(fb_soil)
## Soil depth (mm): 4500
##
## Layer 1 [ 0 to 300 mm ]
## clay (%): 39 silt (%): 29 sand (%): 26 organic matter (%): 6 [ Clay loam ]
## Rock fragment content (%): 50 Macroporosity (%): 7
## Theta WP (%): 25 Theta FC (%): 39 Theta SAT (%): 54 Theta current (%) 39
## Vol. WP (mm): 37 Vol. FC (mm): 58 Vol. SAT (mm): 81 Vol. current (mm): 58
## Temperature (Celsius): NA
##
## Layer 2 [ 300 to 1000 mm ]
## clay (%): 39 silt (%): 32 sand (%): 26 organic matter (%): 3 [ Clay loam ]
## Rock fragment content (%): 65 Macroporosity (%): 7
## Theta WP (%): 24 Theta FC (%): 38 Theta SAT (%): 49 Theta current (%) 38
## Vol. WP (mm): 59 Vol. FC (mm): 93 Vol. SAT (mm): 121 Vol. current (mm): 93
## Temperature (Celsius): NA
##
## Layer 3 [ 1000 to 2000 mm ]
## clay (%): 39 silt (%): 34 sand (%): 26 organic matter (%): 1 [ Clay loam ]
## Rock fragment content (%): 90 Macroporosity (%): 7
## Theta WP (%): 24 Theta FC (%): 37 Theta SAT (%): 47 Theta current (%) 37
## Vol. WP (mm): 24 Vol. FC (mm): 37 Vol. SAT (mm): 47 Vol. current (mm): 37
## Temperature (Celsius): NA
##
## Layer 4 [ 2000 to 4500 mm ]
## clay (%): 39 silt (%): 34 sand (%): 26 organic matter (%): 1 [ Clay loam ]
## Rock fragment content (%): 95 Macroporosity (%): 7
## Theta WP (%): 24 Theta FC (%): 37 Theta SAT (%): 47 Theta current (%) 37
## Vol. WP (mm): 29 Vol. FC (mm): 47 Vol. SAT (mm): 58 Vol. current (mm): 47
## Temperature (Celsius): NA
##
## Total soil saturated capacity (mm): 306
## Total soil water holding capacity (mm): 235
## Total soil extractable water (mm): 106
## Total soil current Volume (mm): 235
## Saturated water depth (mm): NA
We will normally take SpParamsMED as starting point for
species parameters:
data("SpParamsMED")
However, sometimes one may wish to override species defaults with custom values. In the case of FontBlanche there is a table of preferred values for some parameters:
fb$customParams
## Species VCleaf_P12 VCleaf_P50 VCleaf_P88 VCleaf_slope VCstem_P12
## 1 Phillyrea latifolia NA NA NA NA -1.971750
## 2 Pinus halepensis NA NA NA NA -3.707158
## 3 Quercus ilex NA NA NA NA -4.739642
## VCstem_P50 VCstem_P88 VCstem_slope VCroot_P12 VCroot_P50 VCroot_P88
## 1 -6.50 -11.028250 11 NA NA NA
## 2 -4.79 -5.872842 46 -1 -1.741565 -2.301482
## 3 -6.40 -8.060358 30 NA NA NA
## VCroot_slope VCleaf_kmax LeafEPS LeafPI0 LeafAF StemEPS StemPI0 StemAF Gswmin
## 1 NA 3.00 12.38 -2.13 0.5 12.38 -2.13 0.4 0.002
## 2 NA 4.00 5.31 -1.50 0.6 5.00 -1.65 0.4 0.001
## 3 NA 2.63 15.00 -2.50 0.4 15.00 -2.50 0.4 0.002
## Gswmax Gs_P50 Gs_slope Al2As
## 1 0.2200 -2.207094 89.41176 NA
## 2 0.2175 -1.871216 97.43590 631.000
## 3 0.2200 -2.114188 44.70588 1540.671
We can use function modifySpParams() to replace the
values of parameters for the desired traits, leaving the rest
unaltered:
fb_SpParams <- modifySpParams(SpParamsMED, fb$customParams)
Since we are about to run a basic water balance simulation, we
initialize a simulation control parameter list with
transpirationMode = "Granier", i.e.:
fb_control <- defaultControl("Granier")
and we assemble our inputs into a spwbInput object,
using:
fb_x1 <- spwbInput(fb_forest, fb_soil, fb_SpParams, fb_control)
The daily weather data comprises one year:
fb_meteo <- fb$meteoData
nrow(fb_meteo)
## [1] 365
Now, we are ready to launch the simulation:
fb_basic <- spwb(fb_x1, fb_meteo, elevation = 420, latitude = 43.24083)
## Package 'meteoland' [ver. 2.2.2]
## Initial plant water content (mm): 17.1947
## Initial soil water content (mm): 213.886
## Initial snowpack content (mm): 0
## Performing daily simulations
##
## [Year 2014]:............
##
## Final plant water content (mm): 17.1952
## Final soil water content (mm): 236.416
## Final snowpack content (mm): 0
## Change in plant water content (mm): 0.000480745
## Plant water balance result (mm): 0.000480745
## Change in soil water content (mm): 22.5304
## Soil water balance result (mm): 22.5304
## Change in snowpack water content (mm): 0
## Snowpack water balance result (mm): 0
## Water balance components:
## Precipitation (mm) 1066 Rain (mm) 1066 Snow (mm) 0
## Interception (mm) 141 Net rainfall (mm) 925
## Infiltration (mm) 827 Infiltration excess (mm) 98 Saturation excess (mm) 268 Capillarity rise (mm) 0
## Soil evaporation (mm) 26 Herbaceous transpiration (mm) 0 Woody plant transpiration (mm) 313
## Plant extraction from soil (mm) 313 Plant water balance (mm) 0 Hydraulic redistribution (mm) 7
## Runoff (mm) 365 Deep drainage (mm) 198
Surface run-off occurs the same day as precipitation events, whereas deep drainage can last for some days after the event:
g1<-plot(fb_basic)+scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(legend.position = "none")
g2<-plot(fb_basic, "Export")+scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(legend.position = c(0.35,0.60))
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_grid(g1,g2, ncol=1, rel_heights = c(0.5,1))
Precipitation events also generate flows of intercepted water the same day of the event. Evaporation from the bare soil can proceed some days after the event. Transpiration flow is the dominant one in most days, decreasing in summer due to drought.
g1<-plot(fb_basic)+scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(legend.position = "none")
g2<-plot(fb_basic, "Evapotranspiration")+scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(legend.position = c(0.13,0.73))
plot_grid(g1,g2, ncol=1, rel_heights = c(0.5,1))
We can display the dynamics of water potential in different soil layers using:
plot(fb_basic, "SoilPsi")
Tip: Normally, we should expect lower layers to have a less dynamic behaviour, but strange results can occur if, for instance, a large proportion of roots is in deeper layers.
Observations are in element measuredData of the
list:
fb_observed <- fb$measuredData
We can compare the observed vs modelled total evapotranspiration by plotting the two time series:
evaluation_plot(fb_basic, fb_observed, type = "ETR", plotType="dynamics")
It is easy to see that in rainy days the predicted evapotranspiration is much higher than that of the observed data.
We repeat the comparison but excluding the intercepted water from modeled results:
evaluation_plot(fb_basic, fb_observed,
type = "SE+TR", plotType="dynamics")+
theme(legend.position = c(0.8,0.85))
evaluation_plot(fb_basic, fb_observed,
type = "SE+TR", plotType="scatter")
Where we see a reasonably good relationship, but the model tends to underestimate total evapotranspiration during seasons with low evaporative demand.
Function evaluation_stats() allows us to generate
evaluation statistics:
evaluation_stats(fb_basic, fb_observed, type = "SE+TR")
## n Bias Bias.rel MAE MAE.rel r
## 365.00000000 -0.40453685 -30.36201568 0.48384513 36.31440058 0.70498484
## NSE NSE.abs
## 0.03907060 0.02826233
We can now compare the soil moisture content dynamics using:
evaluation_plot(fb_basic, fb_observed, type = "SWC", plotType="dynamics")
The two series have similar shape but not absolute values. This may be an indication that the parameters of the soil water retention curve do not match the data produced by the moisture sensor.
We can calculate evaluation statistics using:
evaluation_stats(fb_basic, fb_observed, type = "SWC")
## n Bias Bias.rel MAE MAE.rel r
## 364.00000000 -0.12212079 -28.47040214 0.12220446 28.48990865 0.91255597
## NSE NSE.abs
## -0.01186163 0.04257990
Note that, despite the bias, the correlation is rather high.
Since we are about to run a advanced water balance simulation, we
initialize a simulation control parameter list with
transpirationMode = "Sperry", i.e.:
fb_control <- defaultControl("Sperry")
and assemble our inputs into a spwbInput object,
using:
fb_x2 <- spwbInput(fb_forest, fb_soil, fb_SpParams, fb_control)
Finally, we launch the simulation (takes ~ 8 seconds):
fb_adv <- spwb(fb_x2, fb_meteo, elevation = 420, latitude = 43.24083)
## Initial plant water content (mm): 31.8864
## Initial soil water content (mm): 213.886
## Initial snowpack content (mm): 0
## Performing daily simulations
##
## [Year 2014]:............
##
## Final plant water content (mm): 31.4636
## Final soil water content (mm): 235.068
## Final snowpack content (mm): 0
## Change in plant water content (mm): -0.422751
## Plant water balance result (mm): -1.97044e-15
## Change in soil water content (mm): 21.1829
## Soil water balance result (mm): 21.1829
## Change in snowpack water content (mm): 0
## Snowpack water balance result (mm): 0
## Water balance components:
## Precipitation (mm) 1066 Rain (mm) 1066 Snow (mm) 0
## Interception (mm) 141 Net rainfall (mm) 925
## Infiltration (mm) 833 Infiltration excess (mm) 92 Saturation excess (mm) 272 Capillarity rise (mm) 0
## Soil evaporation (mm) 21 Herbaceous transpiration (mm) 0 Woody plant transpiration (mm) 323
## Plant extraction from soil (mm) 323 Plant water balance (mm) -0 Hydraulic redistribution (mm) 36
## Runoff (mm) 364 Deep drainage (mm) 195
To compare the performance of the two models with respect to observed data we can calculate the evaluation statistics for soil moisture:
evaluation_stats(fb_basic, fb_observed, type = "SWC")
## n Bias Bias.rel MAE MAE.rel r
## 364.00000000 -0.12212079 -28.47040214 0.12220446 28.48990865 0.91255597
## NSE NSE.abs
## -0.01186163 0.04257990
evaluation_stats(fb_adv, fb_observed, type = "SWC")
## n Bias Bias.rel MAE MAE.rel r
## 364.00000000 -0.13293067 -30.99054249 0.13293067 30.99054249 0.94982546
## NSE NSE.abs
## -0.09692368 -0.04145536
… and for stand evapotranspiration:
evaluation_stats(fb_basic, fb_observed, type = "SE+TR")
## n Bias Bias.rel MAE MAE.rel r
## 365.00000000 -0.40453685 -30.36201568 0.48384513 36.31440058 0.70498484
## NSE NSE.abs
## 0.03907060 0.02826233
evaluation_stats(fb_adv, fb_observed, type = "SE+TR")
## n Bias Bias.rel MAE MAE.rel
## 365.000000000 -0.388630403 -29.168176616 0.470984668 35.349174644
## r NSE NSE.abs
## 0.691737785 -0.009257942 0.054090829
We can compare soil layer moisture dynamics by drawing soil water potentials:
g1<-plot(fb_basic, "SoilPsi", ylim= c(-4,0))+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+theme(legend.position = "none")
g2<-plot(fb_adv, "SoilPsi", ylim= c(-4,0))+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+theme(legend.position = c(0.1,0.47))
plot_grid(g1,g2, ncol=1)
The basic model dries the soil more than the advanced model, which produces a stronger coupling between soil layers because of hydraulic redistribution.
The following shows the daily root water uptake (or release) from different soil layers, and the daily amount of water entering soil layers due to hydraulic redistribution:
g1<-plot(fb_adv, "PlantExtraction")+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+theme(legend.position = "none")
g2<-plot(fb_adv, "HydraulicRedistribution")+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+theme(legend.position = c(0.08,0.5))
plot_grid(g1, g2, rel_heights = c(0.8,0.8), ncol=1)
We can display the transpiration per leaf area unit basis using
"TranspirationPerLeaf".
g1<-plot(fb_basic, "TranspirationPerLeaf", bySpecies = TRUE, ylim = c(0,1.7))+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+theme(legend.position = "none")
g2<-plot(fb_adv, "TranspirationPerLeaf", bySpecies = TRUE, ylim = c(0,1.7))+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+theme(legend.position = c(0.1,0.7))
plot_grid(g1,g2, ncol=1)
The basic model produces higher transpiration than the advanced model.
The following displays the observed and predicted transpiration for Pinus halepensis …
evaluation_plot(fb_adv, fb_observed, cohort = "T2_148", type="E", plotType = "dynamics")+
theme(legend.position = c(0.85,0.83))
evaluation_stats(fb_adv, fb_observed, cohort = "T2_148", type="E")
## n Bias Bias.rel MAE MAE.rel r
## 300.0000000 0.3116129 151.5098933 0.3277188 159.3407793 0.8572620
## NSE NSE.abs
## -11.6137565 -2.3117672
… and Quercus ilex:
evaluation_plot(fb_adv, fb_observed, cohort = "T3_168", type="E", plotType = "dynamics")+
theme(legend.position = c(0.85,0.83))
evaluation_stats(fb_adv, fb_observed, cohort = "T3_168", type="E")
## n Bias Bias.rel MAE MAE.rel
## 309.000000000 -0.005350873 -1.848579701 0.066790546 23.074298510
## r NSE NSE.abs
## 0.888794621 0.764822528 0.544970055
The following plots illustrate that the model simulates a tighter stomatal control for Pinus halepensis.
g1<-plot(fb_adv)+scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(legend.position = "none")
g2<-plot(fb_adv, "LeafPsiRange", bySpecies = TRUE)+
scale_x_date(date_breaks = "1 month", date_labels = "%m")+
theme(legend.position = c(0.1,0.25)) + ylab("Leaf water potential (MPa)")
plot_grid(g1, g2, ncol=1, rel_heights = c(0.4,0.8))
If we compare leaf water potentials against observations
(type = "WP") we obtain a rather good performance for
Q. ilex, but midday water potentials are less well approximated
for P. halepensis.
g1<-evaluation_plot(fb_adv, fb_observed, cohort = "T2_148", type="WP", plotType = "dynamics")
g2<-evaluation_plot(fb_adv, fb_observed, cohort = "T3_168", type="WP", plotType = "dynamics")
plot_grid(g1, g2, ncol=1)
g1<-plot(fb_basic, "PlantStress", bySpecies = TRUE)+
theme(legend.position = "none")
g2<-plot(fb_basic, "StemPLC", bySpecies = TRUE)+
theme(legend.position = c(0.15,0.45))
plot_grid(g1, g2, ncol=1)
g3<-plot(fb_adv, "PlantStress", bySpecies = TRUE)+
theme(legend.position = "none")
g4<-plot(fb_adv, "StemPLC", bySpecies = TRUE)+
theme(legend.position = c(0.15,0.45))
plot_grid(g3, g4, ncol=1)
The basic model seems to overestimate PLC for Pinus halepensis, compared to the advanced model.
This could arise from a difference in the parameters determining PLC or differences in the water potential simulated by both models. We examine the first option using:
plot(fb_basic, "PlantPsi", bySpecies = TRUE)+
theme(legend.position = c(0.15,0.45))
plot(fb_adv, "StemPsi", bySpecies = TRUE)+
theme(legend.position = c(0.15,0.45))
The two models produce similar overall patterns, but the advanced model predicts less negative water potentials for P. halepensis in summer.